
function HL = HessL(B,lambda,X,KD,A,n,W,nu,wnu)

N0 = n(1); N1 = n(2); N2 = n(3); N = N0 + N1 + N2;
[NN,K] = size(X);
R = size(nu,1);

HL = blkdiag(zeros(K+2+NN),2*eye(size(W,1))/W); % Hessian of MPEC Lagrangian

s = X * B(1:K,1) + B(K+1,1) * X(:,K-1) * nu(:,1)' + ...
    B(K+2,1) * X(:,K) * nu(:,2)' + B(K+2+(1:NN),1);
s = exp(s);
s = s ./ (1 + A * s);

xx = zeros(K+2,R,5);
for t = 1:N0
    tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
    ss = reshape(s(tt,:)',1,R,5);
    xx(1:K,:,:) = repmat(reshape(X(tt,:)',K,1,5),1,R,1);
    xx(K+1,:,:) = nu(:,1) * X(tt,K-1)';
    xx(K+2,:,:) = nu(:,2) * X(tt,K)';
    for j = 1:5
        A = zeros(size(HL));
        for k = 1:K+2
            sx = sum(ss.*xx(k:K+2,:,:),3);
            A(k:K+2,k) = (ss(:,:,j) .* (- sum(ss .* xx(k,:,:) .* (xx(k:K+2,:,:) - sx),3) + ...
                (xx(k,:,j) - sx(1,:)) .* (xx(k:K+2,:,j) - sx))) * wnu;
            A(k,k+1:K+2) = A(k+1:K+2,k)';
            A(K+2+tt,k) = A(1+[0,KD+3,2*(KD+3),2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)],k);
            A(k,K+2+tt) = A(K+2+tt,k)';
        end
        A(K+2+tt,K+2+tt) = A(1+[0,KD+3,2*(KD+3),2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)],1+[0,KD+3,2*(KD+3),2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)]);
        HL = HL + lambda.eqnonlin(tt(j)) * A;
    end
end
xx = zeros(K+2,R,4);
for t = 1:N1
    tt = [N0+t,N+N0+t,2*N+N0+N2+N0+t,3*N+N0+N0+t];
    ss = reshape(s(tt,:)',1,R,4);
    xx(1:K,:,:) = repmat(reshape(X(tt,:)',K,1,4),1,R,1);
    xx(K+1,:,:) = nu(:,1) * X(tt,K-1)';
    xx(K+2,:,:) = nu(:,2) * X(tt,K)';
    for j = 1:4
        A = zeros(size(HL));
        for k = 1:K+2
            sx = sum(ss.*xx(k:K+2,:,:),3);
            A(k:K+2,k) = (ss(:,:,j) .* (- sum(ss .* xx(k,:,:) .* (xx(k:K+2,:,:) - sx),3) + ...
                (xx(k,:,j) - sx(1,:)) .* (xx(k:K+2,:,j) - sx))) * wnu;
            A(k,k+1:K+2) = A(k+1:K+2,k)';
            A(K+2+tt,k) = A(2+[0,KD+3,2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)],k);
            A(k,K+2+tt) = A(K+2+tt,k)';
        end
        A(K+2+tt,K+2+tt) = A(2+[0,KD+3,2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)],2+[0,KD+3,2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)]);
        HL = HL + lambda.eqnonlin(tt(j)) * A;
    end
end
for t = 1:N2
    tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,3*N+N0+N0+N1+t];
    ss = reshape(s(tt,:)',1,R,4);
    xx(1:K,:,:) = repmat(reshape(X(tt,:)',K,1,4),1,R,1);
    xx(K+1,:,:) = nu(:,1) * X(tt,K-1)';
    xx(K+2,:,:) = nu(:,2) * X(tt,K)';
    for j = 1:4
        A = zeros(size(HL));
        for k = 1:(K+2)
            sx = sum(ss.*xx(k:K+2,:,:),3);
            A(k:K+2,k) = (ss(:,:,j) .* (- sum(ss .* xx(k,:,:) .* (xx(k:K+2,:,:) - sx),3) + ...
                (xx(k,:,j) - sx(1,:)) .* (xx(k:K+2,:,j) - sx))) * wnu;
            A(k,k+1:K+2) = A(k+1:K+2,k)';
            A(K+2+tt,k) = A(3+[0,KD+3,2*(KD+3)-1,2*(KD+3)+2*(KD+2)],k);
            A(k,K+2+tt) = A(K+2+tt,k)';
        end
        A(K+2+tt,K+2+tt) = A(3+[0,KD+3,2*(KD+3)-1,2*(KD+3)+2*(KD+2)],3+[0,KD+3,2*(KD+3)-1,2*(KD+3)+2*(KD+2)]);
        HL = HL + lambda.eqnonlin(tt(j)) * A;
    end
end

HL = sparse(HL);



